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Abstract 

We present an analytical and numerical study of the velocity distribution function of 
self gravitating collisionless particles, which include dark matter and star clusters. 
We show that the velocity distribution derived through the Eddington's formula 
is identical to the analytical one derived directly from the generalized entropy of 
non-extensive statistical mechanics. This implies that self gravitating collisionless 
structures are to be described by non-extensive thermo-statistics. We identify a 
connection between the density slope of dark matter structures, 7, from p ~ r~'^, and 
the entropic index, q, from the generalized entropy, Sq. Our numerical result confirms 
the analytical findings of earlier studies and clarifies which is the correct connection 
between the density slope and the entropic index. We use this result to conclude 
that from a fundamental statistical mechanics point of view the central density slope 
of self gravitating collisionless dark matter structures is not constrained, and even 
cored dark matter structures are allowed with 7 = 0. We find that the outer density 
slope is bounded by 7 = 10/3. 
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1 Introduction 

The existence of dark matter (DM) has been established beyond any doubt, 
and the formation and evolution of cosmological structures through gravita- 
tional attraction is reasonably well understood. The velocity distribution of the 
resulting structures, however, remains unknown. This is in stark contrast to 
the well established velocity distribution of ideal gases, f{v) ~ exp{—mv'^/2kT), 
which is derived from the extensive Boltzmann-Gibbs entropy. Dark matter, 
however, experiences long-range gravitational interaction, and may therefore 
not obey the rule of extensivity. 
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At the same time numerical simulations provide predictions of steep cen- 
tral density cusps with power law slopes, p ~ r"''', with 7 from 1 to 1.5 
within a few percent of the virial radius of the halo (Navarro, Frenk & White 
1996, Moore et al. 1998). However, recent careful studies (Diemand, Moore 
& Stadel 2004, Reed et al. 2003, Stoehr 2004, Navarro et al. 2004) indicate 
that the resolved region has still not converged on a central density slope. 
Analytical analyses of the Jeans equations seem to indicate that the most 
shallow allowed slope is 7 = 1 (Hansen 2004). Such steep inner numerically 
resolved slopes are, however, not supported by observations. By measuring 
the rotation curve of a galaxy one can in principle determine the density pro- 
file of its DM halo. Low surface brightness galaxies and spirals, where the 
observed dynamics should be DM dominated, seem to show slowly rising ro- 
tation curves (Rubin et al. 1985, Courteau 1997, Palunas & Williams 2000, de 
Blok et al. 2001, de Blok, Bosma & McGaugh 2003, Salucci 2001, Swaters et 
al. 2002, CorbeUi 2003) indicating that these DM halos have constant den- 
sity cores. Galaxy clusters, where baryons can play even less of a role, may 
show a similar discrepancy. Arcs (Sand, Treu & Ellis 2002) and strong lensing 
fits of multiple image configurations and brightnesses (Tyson, Kochanski & 
deir Antonio 1998) also indicate shallow cores in clusters. 

Facing this apparent disagreement between observations and numerical sim- 
ulations it is very important to understand if the pure dark matter central 
density slopes are constrained from a fundamental statistical mechanics point 
of view. We will here combine an analytical and a numerical approach. First we 
show that the dark matter distribution functions are exactly the ones derived 
from non-extensive statistical mechanics. We then find a connection between 
the density slope, 7, and the entropic index, q. A given theoretical bound on 
q will then lead to a bound on 7. We show that the bounds we can imagine 
still allow any density slope, including a core with 7 = 0. The outer slope is 
bounded by 10/3, which is close to the findings of N-body codes. Our find- 
ings can be very useful and directly applicable if one can envisage stronger 
theoretical bounds on q. 



2 Eddington's formula 

Let us first derive the actual velocity distribution function of self gravitating 
coUisionless structures numerically. 

We will be considering the simplest possible dark matter structures (actu- 
ally any coUisionless system, including star clusters), which are spherical and 
isotropic. According to recent numerical N-body results this is a reasonably 
good approximation in the central equilibrated part of dark matter structures. 
We will argue later that statistical mechanics in a natural way forces a system 



2 



towards isotropy. Any given density profile, p(r), can be inverted (Eddington 
1916) to give the particle velocity distribution function, f{E), through 



where ^'(r) is the relative potential as a function of radius, and E is the 
relative energy, E = — mv^/2. For recent applications and technical de- 
tails see (Binney 1982, Kazantzidis et al. 2004). For the numerical inversion 
we typically consider structures with 10 orders of magnitude in radius, and 
always confirm that the resulting velocity distributions have converged. We 
consider only isotropic structures and leave non-isotropy for a future analysis. 
For simplicity we construct very simple structures which depend only on one 
free parameter, namely the density slope 7, from p(r) ~ r~^. The resulting 
velocity distribution functions are presented as (coloured) symbols in figure 1 
for the slopes 7 = 1.0, 1.5, 2 and 2.3. 

It is important to note that for 7 > 2 there is a maximal velocity. This is a real 
feature, purely related to the density profile, and is not simply the result of 
a finite structure. This shows directly that the velocity distribution function 
cannot be an exponential or any sum of exponentials, since exponential func- 
tions have tails going to infinite energy. We have thus shown explicitely that 
the often discussed sum of exponentials (Lynden-Bell 1967, KuU, Treumann 
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Fig. 1. The velocity distribution functions. The (coloured) symbols are from invert- 
ing a given density distribution according to Eddington's formula, for the density 
slopes of 7 = 1.0, 1.5, 2.0 and 2.3. The solid lines are the theoretical formulae from 
Tsallis statistics, using entropic indices of g = 1.39, 1.21, 1.0 and 0.86. The curves are 
arbitrarily normalized, and use different velocity dispersion, a^, in order to make the 
figure more readable and the difference in shapes more evident. The labels describe 
the lines from left to right. 
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& Bohringer 1997, Iguchi et al. 2004) cannot be correct, and instead a more 
general functional form for the velocity distribution function must be sought. 
We will now turn to non-extensive thermo-statistics and show that the correct 
functional form for the velocity distribution function is easy to derive. 



3 Tsallis statistics 



Statistical mechanics for classical gases can be derived from the Boltzmann- 
Gibbs assumption for the entropy, Sbg — ~k^Pi • Inpj, where Pi is the prob- 
ability for a given particle to be in the state i, and the sum is over all states. 
For normal gases the probability, p{v), coincides with the velocity distribu- 
tion function, f{v). This classical statistics can be generalized to Tsallis (or 
non-extensive) statistics (TsaUis 1988), which depends on the entropic index 
Q 

Sq = -k^Pi -InqPi, (2) 



where the q- logarithm is defined by, lugp — — 1)/(1 — g) , and ior q — 1 the 
normal Boltzmann-Gibbs entropy is recovered, Sbg = •S'l- The probabilities 
still obey, J2Pi = Ij while the particle distribution function is now given by 
f{v) = p^{v). Thus, for g < 1 one privileges rare events, whereas q > 1 
privileges common events. One also sees directly that entropy is maximized 
for isotropic distributions, which supports our original isotropy assumption. 
For a summary of applications see (Tsallis 1999), and for up to date list of 
references see http://tsallis.cat.cbpf.br/biblio.htm. 

Average values are calculated through the particle distribution function, and 
one e.g. has the mean energy (Tsallis, Mendes & Plastino 1998) 

U, = ^, (3) 

Cp 



where Cp = J2Pi, and Ei are the energy eigenvalues. Optimization of the 
entropy in eq. (2) under the constraints leads to the probability (Silva, Plastino 
& Lima 1998, TsaUis, Prato & Plastino 2003) 

_ [l-{l-q)f3,{E,-U,)f^'-'^^ 



where Zg normalizes the probabilities, (3q = (3/cp, and (3 is the optimization 
Lagrange multiplier associated with the average energy. Adding a constant 
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energy, eo, to all the energy eigenvalues leads to Uq Uq + eo, which leaves 
all the probabilities, Pi, invariant (Tsallis, Mendes & Plastino 1998). When a 
is positive, where 



a = 1 + (1 - q)PqUq , 



(5) 



then eq. (4) can be written as 



{l-{l-q){PJa)E,) 



1/(1-9) 



Pi = 



(6) 



On figure 1 we plot 4 solid (black) lines, corresponding to q — 1.39, 1.21, 1.0 
and 0.86, and well chosen /3's. 

It is clear from figure 1 that the velocity distribution function, f{v), inverted 
from a density distribution function with slope 7 = 2.3, where p ~ r~'^, is 
extremely well fitted with the non-extensive form in eq. (6) using q = 0.86. 
Similarly, the Maxwell distribution (the exponential, g = 1) is recovered with 
an isothermal sphere (with 7 = 2), and finally, the j — 1.5 is well fitted with 
q — 1.21. One should keep in mind that the probabilities, p{v), are related to 
the physical distributions through f{v) — p^{v). 

We consider this excellent agreement in figure 1 as a strong indicator that equi- 
librium self gravitating structures indeed are to be described by non-extensive 
statistical mechanics. We are thus giving theoretical support to the use of Tsal- 
lis entropy (and the resulting velocity distribution function) when analysing 
gravitational structures. Such analysis was first conducted in (Lavagno et 
al. 1998), where an entropic index of g ~ 0.23 was found to describe accu- 
rately the peculiar velocity of galaxy clusters. It is worth emphasizing that 
the maximal velocity for < 1, t^max = \j'2a./m(5q{l — q) is an inherent prop- 
erty of the non-extensive thermo-statistics (for q <1) and is not related to any 
finite size of the system. We also note that the a in eq. (5) may be important 
for cluster temperature observations (Hansen 2004). 

The structures we are considering are polytropes, which are described by the 
polytropic index, n, making the connection between pressure and density, p ~ 
pi+i/n^ It was first shown analytically by (Plastino & Plastino 1993) that these 
structures are exactly described by non-extensive statistical mechanics, and a 
connection was found between n and the entropic index, q. Thus, the finding 
of this section is a numerical confirmation of the original work of (Plastino 
& Plastino 1993), which says that the distribution function for polytropes is 
infact given by eq. (6). As we will discuss below, the actual connection between 
n and q is not clear from the literature, and our independent results are useful 
in distinguishing which of the analytical results is correct. 
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Naturally, for a more general density profile where 7 is not constant, the actual 
velocity distribution function will be a sum over terms of the shape given by 
eq. (6), with f{v) — p'^{v). One can only approximate the velocity distribution 
function with one term of the form given in eq. (6) when 7(r) varies sufficiently 
slowly with radius. 



4 A connection between the density slope and the entropic index? 



We can now proceed as in the previous section with various density profiles, 
and find the corresponding q. The result is presented in figure 2. The fitted 
g's are accurate at the few percent level. This figure is presented only in the 
range < 7 < 2.5, since steeper slopes lead to numerical problems in our 
simple inversion of Eddington's formula. This is not a problem since we are 
interested in exactly this range in density profiles. As pointed out in section 
3, our numerical finding that the self-gravitating structures are described by 
non-extensive statistical mechanics was first shown analytically by (Plastino & 
Plastino 1993). There seems, however, to be some confusion in the literature, as 
to what the connection between n and q really is. Polytropes follow the Lane- 
Emden equation, and for power- law density profiles one thus has, n = 7/(7 — 
2). Now, there are 3 different connections between n and q frequently used 
in the literature (see e.g. (Taruya & Sakagami 2003, Chavanis 2003, Tsallis, 
Prato & Plastino 2003) for different expressions) and we can therefore use our 
independent numerical method to distinguish which is the correct one. We 
find that the formula derived in (Tsallis, Prato & Plastino 2003, Lima & de 
Souza 2004, Silva & Alcaniz 2003), n = (5 - 3g)/(2(l - q)) fits our results 
very well. We can therefore conclude that the correct connection between the 
density slope, 7, and the entropic index, q, is given by 

? = n — ■ iv 

D — 7 



This is shown as the thin (green) solid line in figure 2. We note that this equa- 
tion may play an important role in explaining the observed correlation be- 
tween the density slope and the anisotropy of dark matter structures (Hansen 
& Moore 2004). 

It is interesting to ask the question: which are the theoretical constraints on 
the entropic index ql We saw above that eq. (6) is only valid for positive 
a. Actually the expression in eq. (6) is also valid for negative a, however, 
a negative effective temperature (a negative velocity dispersion in our case), 
f^q ~ Pq/'^ is difficult to interpret for cosmological structures. Since we find 
numerically that UqPq = 3/2, one sees that a = for g = 5/3. The same 
conclusion, namely q < 5/3, was reached by (Boghosian 1999) by demanding a 
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positive proportion between temperature and internal energy; and also reached 
by (Silva & Alcaniz 2003) through an analysis of the velocity dispersion. 

We see in figure 2, that g = 5/3 corresponds to a cored dark matter distribution 
with 7 = 0. Thus all dark matter density slopes are allowed from a fundamen- 
tal statistical mechanics point of view. Our findings are directly applicable 
if one can imagine a stronger theoretical bound on q. E.g. (Boghosian 1999) 
considered the ideal g and found that positiveness of the thermal con- 

ductivity leads to q < 7/5. Such a bound would according to our figure 2 
imply a bound on the central density profile of 7 > 1, in agreement with nu- 
merical and analytical suggestions (Diemand, Moore & Stadel 2004, Reed et 
al. 2003, Navarro et al. 2004, Hansen 2004). It has indeed been suggested (Xiao, 
Sun & Hao 2004) that the numerical methods induce an artificial conductivity, 
which would then imply a central slope of —1. 

The lower bound on q is given by g > (Lima, Silva & Plastino 2001), which 
implies a bound on the dark matter profile of 7 < 10/3, when using eq. (7). 
This is very close to the findings of numerical simulations, which is 7 < 3. 
Further, if one would have a bound of g > 1/3 (as suggested by (Lima, Silva 

6 Santos 2001)), then one would have a bound on the outer density profile of 

7 < 3. 
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Fig. 2. A given density profile, p ~ r"'*', gives a specific distribution function, f{v), 
which can be fitted with two free parameters, namely the entropic index, q, and 
the velocity dispersion, a^. We present the resulting g as a function of the density 
slope, 7. The theoretical limit of (7 < 5/3 is presented with a dashed line. The thin 
(green) solid line is the curve given in eq. (7). 
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5 Conclusions 



We compare the velocity distribution function for self gravitating coUisionless 
structures (e.g. dark matter or star clusters) with the analytical expectation 
from non-extensive statistical mechanics, and find excellent agreement be- 
tween the two. This is a strong indication that self gravitating systems of 
coUisionless particles must be described by non-extensive statistical mechan- 
ics. 

We identify a connection between the density power slope, 7, from p ~ r"''', 
and the entropic index, g, from the generalized entropy, Sq. We show that q 
is only bounded from above by q < 5/3, which implies that all dark matter 
central slopes are allowed from this fundamental statistical mechanics point 
of view, including central cores with 7 = 0. We show that the outer density 
slope is bounded by 7 < 10/3. If one can imagine stronger theoretical bounds 
on q then our findings are directly applicable to give stronger bounds on the 
allowed density slope. 
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